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Abstract: We develop a Bayesian nonparametric extension of the popular Plackett-Luce choice model that 
can handle an infinite number of choice items. Our framework is based on the theory of random atomic 
measures, with the prior specified by a gamma process. We derive a posterior characterization and a simple 
and effective Gibbs sampler for posterior simulation. We develop a time-varying extension of our model, 
and apply it to the New York Times lists of weekly bestselling books. 
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Modeles bayesiens non parametriques pour les donnees de rang 



Resume : On s'interesse dans ce rapport a une extension bayesienne non parametrique du modele 
de Plackett-Luce pour les donnees de rang, pouvant traiter un nombre potentiellement infini d'elements. 
Notre cadre se base sur la theorie des mesures completement aleatoires, avec comme a priori un processus 
de gamma. Nous derivons une caracterisation de la loi a posteriori et un echantillonneur de Gibbs simple 
pour approcher la loi a posteriori. Nous developpons egalement une version dynamique de notre modele, 
et l'appliquons aux listes hebdomadaires des 20 meilleures ventes du New York Times. 

Mots-cles : Modeles de choix, modele de Bradley-Terry generalise, modele de Plackett-Luce, processus 
de gamma, methodes de Monte Carlo par chaine de Markov 
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1 Introduction 

Data in the form of partial rankings, i.e. in terms of an ordered list of the top-m items, arise in many 
contexts. For example, in this paper we consider datasets consisting of the top 20 bestselling books 
as published each week by the New York Times. The Plackett-Luce model flT] [2j is a popular model for 
modeling such partial rankings of a finite collection of M items. It has found many applications, including 
choice modeling [3|, sport ranking |4|, and voting Q. (6] Chap. 9] provides detailed discussions on the 
statistical foundations of this model. 

In the Plackett-Luce model, each item k £ [M] — {1, . . . , M} is assigned a positive rating parameter 
Wk, which represents the desirability or rating of a product in the case of choice modeling, or the skill of 
a player in sport rankings. The Plackett-Luce model assumes the following generative story for a top-m 
list p = (pi, . . . , p m ) of items pi £ [M]: At each stage i — 1, . . . , m, an item is chosen to be the ith item 
in the list from among the items that have not yet appeared, with the probability that pi is selected being 
proportional to its desirability w Pi . The overall probability of a given partial ranking p is then: 

m 



W 



Pi, 



with the denominator in ([T} being the sum over all items not yet selected at stage i. 

In many situations the collection of available items can be very large and potentially unknown. In this 
case, a nonparametric approach can be sensible, where the pool of items is assumed to be infinite and the 
model allows for the possibility of items not observed in previous top-m lists to appear in new ones. In 
this paper we propose such a Bayesian nonparametric Plackett-Luce model. Our approach is built upon 
recent work on Bayesian inference for the (finite) Plackett-Luce model and its extensions |7 8 9]. Our 
model assumes the existence of an infinite pool of items {Xkl^Li, each with its own rating parameter, 
{ w k}kLi' The probability of a top-m list of items, say (X pi , . . . , X Pm ), is then a direct extension of the 
finite case (fTh: 

m 

p(x pi ,. . . , x p j = n — — wp \ — (2) 

To formalize the framework, a natural representation to encapsulate the pool of items along with their 
ratings is using an atomic measure: 



G = J2 w kS Xk (3) 



k=X 

Using this representation, note that the top item X Pl in our list is simply a draw from the probability 
measure obtained by normalizing G, while subsequent items in the top-m list are draws from probability 
measures obtained by first removing from G the atoms corresponding to previously picked items and 
normalizing. Described this way, it is clear that the Plackett-Luce model is basically a partial size-biased 
permutation of the atoms in G [10], and the existing machinery of random measures and exchangeable 
random partitions [ 1 1 1 can be brought to bear on our problem. 

In particular, in Section|2]we will use a gamma process as the prior over the atomic measure G. This is 
a completely random measure 1 12 1 with gamma marginals, such that the corresponding normalized prob- 
ability measure is a Dirichlet process. We will show that with the introduction of a suitable set of auxiliary 
variables, we can characterize the posterior law of G given observations of top-m lists distributed accord- 
ing to Q. A simple Gibbs sampler can then be derived to simulate from the posterior distribution. In 
Section|3]we develop a time-varying extension of our model and derive a simple and effective Gibbs sam- 
pler for posterior simulation. In Section [4] we apply our time-varying Bayesian nonparametric Plackett- 
Luce model to the aforementioned New York Times bestsellers datasets, and conclude in Section[5] 
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2 A Bayesian nonparametric model for partial ranking 

We start this section by briefly describing a Bayesian approach to inference in finite Plackett-Luce models 
|9|, and taking the infinite limit to arrive at the nonparametric model. This will give good intuitions for 
how the model operates, before we rederive the same nonparametric model more formally using gamma 
processes. Throughout this paper we will suppose that our data consists of L partial rankings, with 
Pi = (piii ■ ■ ■ , Pim) for t € [L]. For notational simplicity we assume that all the partial rankings are 
length m. 

2.1 Finite Plackett-Luce model with gamma prior 

Suppose we have M choice items, with item k £ [A I] having a positive desirability parameter Wk- A 
partial ranking pi = (pa, . . . , pi m ) can be constructed generatively by picking the ith item pu at the 
ith stage for i — 1, . . . , m, with probability proportional to w PH as in ([!}. An alternative Thurstonian 
interpretation, which will be important in the following, is as follows: For each item k let zgk ~ Exp(u>j,) 
be exponentially distributed with rate Wk- Thinking of z«, as the arrival time of item k in a race, let pu 
be the index of the ith item to arrive (the ith smallest value among (zgk)kLi)- The resulting probability 
of pi can then be shown to still be ([T]). In this interpretation (zgk) can be understood as latent variables, 
and the EM algorithm can be applied to derive an algorithm to find a ML parameter setting for (wk)k=i 
given multiple partial rankings. Unfortunately the posterior distribution of (zgk) given pg is difficult to 
compute directly, so we instead consider an alternative parameterization: Let Zm = z pu — z pe ._ 1 be the 
waiting time for the ith item to arrive after the i — 1th item (with z pta defined to be 0). Then it can be 
shown that the joint probability is: 

L rn 

p(o>t)u, (^=M = iiK)£ii) = nn% ex p (-^ (sf =1 - es «»«,)) (4) 

i=\ i=l 

Note that the posterior of (Za)"^ is simply factorized with Zu\p, w ~ ExpQ3 fe=1 Wk ~ S}=i w ptj)' 
and the ML parameter setting can be easily derived as well. Taking a further step, we note that a factorized 
gamma prior over (wk) is conjugate to Q, say Wk ~ Gamma(-g, r) with hyperparameters a, r > 0. 
Now Bayesian inference can be carried out either with a VB EM algorithm, or a Gibbs sampler. In this 
paper we shall consider only Gibbs sampling algorithms. In this case the parameter updates are of the 
form 

w k \{pi), (Zti), (wk')k'^k ~ Gamma (j^ +n k ,r + J2i=i J2"Li SimZti) (5) 

where rik is the number of occurrences of item k among the observed partial rankings, and 8uk = if 
there is a j < i with pej = k and 1 otherwise. These terms arise by regrouping those in the exponential 
in Q. 

A nonparametric Plackett-Luce model can now be easily derived by taking the limit as the number 
of choice items M —> oo. For those items k that have appeared among the observed partial rankings, 
the limiting conditional distribution <|5j is well defined since rik > 0. For items that did not appear in 
the observations, |5]) becomes degenerate at 0. Instead we can define = Ylk-n =o Wk t0 ^ e tne ioia ^ 
desirability among all infinitely many previously unobserved items, and show that 

W*\(pt), {Zh), (w k )k:n k >o ~ Gamma (a, r + J2i=i YJiLi Z £i) ( 6 ) 

The Gibbs sampler thus alternates between updating (Zn), and updating the ratings of the observed items 
(vJk)k:n k >o and of the unobserved ones w*. This nonparametric model allows us to estimate the proba- 
bility of seeing new items appearing in future partial rankings in a consistent manner. While intuitive, this 
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derivation is ad hoc in the sense that it arises as the infinite limit of the Gibbs sampler for finite models, 
and is unsatisfying as it did not directly capture the structure of the underlying infinite dimensional object, 
which we will show in the next subsection to be a gamma process. 



2.2 A Bayesian nonparametric Plackett-Luce model 

Let X be a measurable space of choice items. A gamma process is a completely random measure over 
X with gamma marginals. Specifically, it is a random atomic measure of the form ([3]), such that for 
each measurable subset A, the (random) mass G(A) is gamma distributed. Assuming that G has no 
fixed atoms (that is, for each element x <E X we have G({x}) = with probability one) and that the 
atom locations {X^} are independent of their masses {wk}, it can be shown that such a random measure 
can be constructed as follows: each is iid according to a base distribution H (which we assume 
is non-atomic with density h(x)), while the set of masses {wk} is distributed according to a Poisson 
process over K + with intensity X(w) = aw~ 1 e~ WT where a > is the concentration parameter and 
r > the inverse scale. We write this as G ~ T(a,T, H). Under this parametrization, we have that 
G(A) ~ G&mm.&{aH{A), r). 

Each atom is a choice item, with its mass > corresponding to the desirability parameter. The 
Thurstonian view described in the finite model can be easily extended to the nonparametric one, where a 
partial ranking (X pei . . . X Pfm ) can be generated as the first m items to arrive in a race. In particular, for 
each atom X^ let Zfk ~ Exp(wfc) be the time of arrival of Xk and X pei the zth item to arrive. The first 
m items to arrive (X pn . . . X ptm ) then constitutes our top-m list, with probability as given in ©. Again 
reparametrizing using inter-arrival durations, let Zn — z pu — z pui for i = 1,2, .. . (with z po = 0). 
Then the joint probability is: 

P((Xp«)iLi, {Zu)? =1 \G) = P((z Pfl . . . z pe J, and z ek > z Pfm for all k $ {p a , Plm \) (7) 




Marginalizing out (Zu)il i gives the probability of {X pti ) r [L 1 in p}. Further, conditional on pi it is seen 
that the inter-arrival durations Zm . . . Zg m are mutually independent and exponentially distributed: 

/ oo i— 1 >. 

Z tl I (X PH )T=i , G ~ Exp ( w k — ^ w pej ) (8) 
^ fe=i j=i ' 

The above construction is depicted on Figure [TJleft). We visualize on right some top-m lists generated 
from the model, with r = 1 and different values of a. 



2.3 Posterior characterization 

Consider a number L of partial rankings, with the it\\ list denoted Yt = (Yn . . . Yg me ) , for £ E [L]. 
While previously our top-m list (X Pl . . . X Pm ) consists of an ordered list of the atoms in G. Here G 
is unobserved and (Ygx ■ ■ ■ Yg me ) is simply a list of observed choice items, which is why they were not 
expressed as an ordered list of atoms in G. The task here is then to characterize the posterior law of G 
under a gamma process prior and supposing that the observed partial rankings were drawn iid from the 
nonparametric Plackett-Luce model given G. Re-expressing the conditional distribution |2]) of Yg given 
G, we have: 
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Figure 1: Bayesian nonparametric Plackett-Luce model. Left: 
G and U = J2 k u k Sx k where u k = — log(z/ c ). The top-3 
ranking is (pi, p2, P3)- Right: Visualization of top-5 rankings 
with rows corresponding to different rankings and columns to 
items sorted by size biased order. A lighter shade corresponds 
to a higher rank. Each figure is for a different G, with a = 
.1,1,3. 



As before, for each I, we will also introduce a set of auxiliary variables Zi = (Zgi . . . Zg me ) (the inter- 
arrival times) that are conditionally mutually independent given G and Yg, with: 

Zti\Yt, G ~ Exp(G(X\{y a , . . . , y«_i})) (10) 

The joint probability of the item lists and auxiliary variables is then (c.f. |7|): 

L mj 

P((Yt, Z g )\ =l \G) = J] II G ({ y «» exp(-Z tt G(X\{Y tl , . . . , l^-i})) (11) 



Note that under the generative process described in Section 2.2 there is positive probability that an item 
appearing in a list Ye appears in another list Yy with t' ^ £. Denote the unique items among all L lists by 
XI . . . Xft, and for each fc = 1, . . . , K let rife be the number of occurrences of X^ among the item lists. 
Finally define occurrence indicators 

fo = (12) 

I 1 otherwise. 

i.e. 8uk is the indicator of the occurence that item X k does not appear at a rank lower than i in the £th 
list. Then the joint probability under the nonparametric Plackett-Luce model is: 

K L mi 

P((Y e ,Z e )f =1 \G) = l[G({X* k }r x HUcM-ZuGiXMYa,...^^})) 

k=l £=li=l 

= exp ( -G(X) Z A II G (i x n) nk exp ( -G({X* k }) ]T(& fc - 1)Z U j 

V U / fc=l V ft / 

(13) 

Taking expectation of ( fT3j ) with respect to G using the Palm formula gives: 

Theorem 1 77ie marginal probability of the L partial rankings and auxiliary variables is: 

P((Y e ,Z e )f =1 ) = e -^£« z «) n h(X* k )Jn k ,Y,5ukZ e ?j (14) 

fe=i ^ a ' 
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where tp(z) is the Laplace transform of X, 

Mz) = -logE [e- zG ^l = f \(w){l-e- zw )dw = a\og(l + -) (15) 

and k(ji, z) is the nth moment of the exponentially tilted Levy intensity X(w)e~ zw : 

K (n,z)= [ \(w)w n e~™dw = ° r(n) (16) 

Details are given in the appendix. Another application of the Palm formula now allows us to derive a 
posterior characterisation of G: 

Theorem 2 Given the observations and associated auxiliary variables (Ye, Ze)\—\> the posterior law of 
G is also a gamma process, but with atoms with both fixed and random locations. Specifically, 

K 

G\(Ye,Z e )f =1 = G*+J2<5 x * (17) 

k=l 

where G* and w*,..., w* K are mutually independent. The law of G* is still a gamma process, 

G*\(Xe,Ze)f =1 ~r(a,T*,h) r*= T + Y J Ze l (18) 

u 

while the masses have distributions, 

w l\( Y ei z t)e=i ~ Gamma f rtfe,r + 5 m Zti ) (19) 



2.4 Gibbs sampling 

Given the results of the previous section, a simple Gibbs sampler can now be derived, where all the 
conditionals are of known analytic form. In particular, we will integrate out all of G* except for its total 
mass w'* = G*(K). This leaves the latent variables to consist of the masses w*, (w%.) and the auxiliary 
variables (Zu)- The update for Ze% is given by ( fTO] ), while those for the masses are given in Theorem|2] 

Gibbs update for ZeC Z tt |rest ~ Exp (w* + J2 k 8eikW* h ) (20) 

Gibbs update for u>£: w^lrest ~ Gamma (rik, r + J2u ^UkZei) (21) 

Gibbs update for w* : w* |rest ~ Gamma (a, t + J2 ei Zu) (22) 

Note that the auxiliary variables are conditionally independent given the masses and vice versa. Hyperpa- 
rameters of the gamma process can be simply derived from the joint distribution in Theorem[T] Since the 
marginal probability of the partial rankings is invariant to rescaling of the masses, it is sufficient to keep 
t fixed at 1. As for a, if a Gamma(a, b) prior is placed on it, its conditional distribution is still gamma: 

Gibbs update for a: a|rest - Gamma (a + K, b + log (l + (23) 

Note that this update was derived with w* marginalized out, so after an update to a it is necessary to 
immediately update w'* via ( |22] i before proceeding to update other variables. 
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3 Dynamic Bayesian nonparametric ranking models 

In this section we develop an extension of the Bayesian nonparametric Plackett-Luce model to model 
time-varying rankings, where the rating parameters of items may change smoothly over time and reflected 
in a changing series of rankings. Given a series of times indexed by t = 1,2,..., we may model the 



rankings at time t using a gamma process distributed random measure Gt as in Section 2.2 with Markov 



dependence among the sequence of measures (Gt) enabling dependence among the rankings over time. 

3.1 Pitt- Walker dependence model 

We will construct a dependent sequence (Gt) which marginally follow a gamma process T(a,r, H) using 
the construction of lfl3ll . Suppose Gt ~ T(a, t, H). Since Gt is atomic, we can write it in the form: 

oo 

Gt = £«*fcfcr tt (24) 
fc=i 

Define a random measure Ct with conditional law: 

oo 

Ct\G t =};CtkSx tk c t k I G t ~ PoissonOi«7 tfc ) (25) 

fc=i 

where 4> t > is a dependence parameter. Using the same method as in Section |2~3] we can show: 
Proposition 3 Suppose the law of Gt is T(a, t, H). The conditional law of Gt given Ct is then: 

oo 

G t = G* t +J2 w *tkSx tk (26) 
fc=i 

where G* t and (if^O^Li are a ^ mutually independent. The law ofG* t is given by a gamma process, while 
the masses are conditionally gamma, 

G* t \C t ~T(a,T + ^ t ,H) u4|C t ~Garnma(ct fc ,T + t ) (27) 

The idea of lfl3l is to define the conditional law of G t +\ given G t and C t to coincide with the 
conditional law of Gt given Ct as in Proposition^ In other words, define 

oo 

G t +\ = G* t+l + ^ w t+i,k$x tk (28) 

k=l 

where G* +1 ~ T(a, t + <fit,H) and w t +i,k ~ Gamma(c f t,T + <fi t ) are mutually independent. If the 
prior law of G t is T(a, r, H), the marginal law of G t +i will be T(a, t, H) as well when both G t and 
C t are marginalized out, thus maintaining a form of stationarity. Further, although we have described the 
process in order of increasing t, the joint law of Gt, Ct, G*+i can equivalently be described in the reverse 
order with the same conditional laws as above. Note that if c t k — 0, the conditional distribution of w t +i.k 
will be degenerate at 0. Hence Gt+i has an atom at X t k if and only if Ct has an atom at X t k, that is, 
if Ctk > 0. In addition, it also has atoms (those in G*,{) where Ct does not (nor does Gt). Finally, the 
parameter cp t can be interpreted as controlling the strength of dependence between G t +\ and G t - Indeed 
it can be shown that 

E[G i+1 \G t ] - -^—Gt + -^—H. (29) 
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Another measure of dependence can be gleaned by examining the "lifetime" of an atom. Suppose X 
is an atom in G\ with mass w > 0. The probability that X is an atom in C2 with positive mass is 
1 — exp(— (j>iw), in which case it has positive mass in G2 as well. Conversely, once it is not an atom, 
it will never be an atom in the future since the base distribution H is non-atomic. The lifetime of the 
atom is then the smallest t such that it is no longer an atom. We can show by induction that: (details in 
appendix) 

Proposition 4 The probability that an atom X in G\ with mass w > is dead at time t is given by 

P(G t ({X}) = 0\w)=en>(-yt\iw) 
where y t \i can be obtained by the recurrence Dt\t-1 — (f>t— 1 an d Vt\s—\ — $ Vt [ 3 ^'r+y { • 

3.2 Posterior characterization and Gibbs sampling 

Assume for simplicity that at each time step t = 1, . . . , T we observe one top-m list Y t — (Y tl , . . . , Y tm ) 
(it trivially extends to multiple partial rankings of differing sizes). We extend the results of the previous 
section in characterizing the posterior and developing a Gibbs sampler for the dynamical model. 

Since each observed item at time t has to be an atom in its corresponding random measure Gt, and 
atoms in Gt can propagate to neighboring random measures via the Pitt- Walker dependence model, we 
conclude that the set of all observed items (through all times) has to include all fixed atoms in the posterior 
of Gt- Thus let X* = (X%), k = 1, . . . ,K be the set of unique items observed in Y\, . . . , Yt, let 
n t k € {0, 1} be the number of times the item X£ appears at time t, and let p t be defined as Y t — 
(X* , . . • , X* m ). We write the masses of the fixed atoms as wtk — Gt({X£}), while the total mass of 
all other random atoms is denoted wt* = Gt(X\X*). Note that wtk has to be positive on a random 
contiguous interval of time that includes all observations of X£ — it's lifetime — but is zero outside of the 
interval. We also write c tk — C t {{X%,}) and c t * = C t (X\X*). As before, we introduce, for t = 1, . . . , T 
and i — 1, . . . , m, latent variables 

( K ^ \ 

Z tl ~ Exp ( w t * + X! Wtk ~ Wt Pi ) ( 3 °) 
^ k=i 3=1 ' 

Each iteration of the Gibbs sampler then proceeds as follows (details in appendix). The latent variables 
iZ t i) are updated as above. Conditioned on the latent variables (Z t i), (ctk) and (q*), we update the 
masses (w tk ), which are independent and gamma distributed since all likelihoods are of gamma form. 
Note that the total masses (Gt(X)) are not likelihood identifiable, so we introduce an extra step to improve 
mixing by sampling them from the prior (integrating out (ctk), (ct*)), scaling all masses along with it. 
Directly after this step we update (c t k), (ct* ) ■ We update a along with the random masses (w t * ) and (c t * ) 
efficiently using a forward-backward recursion. Finally, the dependence parameters (<fr t ) are updated. 



3.3 Continuous time formulation using superprocesses 

The dynamic model described in the previous section is formulated for discrete time data. When the 
time interval between ranking observations is not constant, it is desirable to work with dynamic models 
evolving over continuous-time instead, with the underlying random measures (Gt) defined over all (el, 
but with observations at a discrete set of times t\ < t<i < • • • . Here we propose a continuous-time model 
based on the Dawson-Watanabe superprocess lfl4l[T5ll (see also lfl6l[T7l[T8l[T9l ). This is a diffusion on 
the space of measures with the gamma process T(a, r, H ) as its equilibrium distribution. It is defined by 
a generator 

c = * (/ G ^mxy +a J H ^dix) - T S G ^mx)) 
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Figure 2: Sample path drawn from the Dawson-Watanabe superprocess. Each colour represents an atom, 
with height being its (varying) mass. Left shows (Gt) and right (Gt/Gt(X)), a Fleming- Viot process. 

with £ parametrizing the rate of evolution. Figure [2] gives a sample path, where we see that it is con- 
tinuous but non-differentiable. For efficient inference, it is desirable to be able to integrate out all Gf's 
except those G tl ,Gt 2 , ... at observation times. An advantage to using the Dawson-Watanabe super- 
process is that, the conditional distribution of Gt s given Gt s _ 1 is remarkably simple ll20l . In particu- 
lar it is simply given by the discrete-time process of the previous section with dependence parameter 
0tji„_i — t;(* s -*s_i) — ■ Thus the inference algorithm developed previously is directly applicable to 
the continuous-time model too. 

4 Experiments 

We apply the discrete-time dynamic Plackett-Luce model to the New York Times bestsellers data. These 
consist of the weekly top-20 best-sellers list from June 2008 to April 2012 in various categories. We 
consider here the categories paperback nonfiction (PN) and hardcover fiction (HF), for which respectively 
249 and 916 books appear at least once in the top-20 lists over the 200 weeks. We consider that the 
correlation parameter <j>t = <fi is constant over time, and assign flat improper priors p(a) oc 1/a and 
p{4>) oc 1/0. In order to take into account the publication date of a book, we do not consider books in the 
likelihood before their first appearance in a list. We run the Gibbs sampler with 10000 burn-in iterations 
followed by 10000 samples. Mean normalized weights for the more popular books in both categories are 
shown in Figure [5] 

The model is able to estimate the weights associated to each book that appeared at least once, as well 
as the total weight associated to all other books, i.e. the probability that a new book enters at the first rank 
in the list, represented by the black curve. Moreover, the Bayesian approach enables us to have a measure 
of the uncertainty on the weights. The hardcover fiction category is characterized by rapid changes in 
successive lists, compared to the paperback nonfiction. This is quantified by the estimated value of the 
parameter 4>, which are respectively 85 ± 20 and 140 ± 40 for PN and HF. The estimated values of the 
shape parameter a are 7 ± 1.5 and 2 ± 1 respectively. 

5 Discussion 

We have proposed a Bayesian nonparametric Plackett-Luce model for ranked data. Our approach is based 
on the theory of atomic random measures, where we showed that the Plackett-Luce generative model 
corresponds exactly to a size-biased permutation of the atoms in the random measure. We characterized 
the posterior distribution, and derived a simple MCMC sampling algorithm for posterior simulation. 
Our approach can be see as a multi-stage generalization of posterior inference in normalized random 
measures ll2Tll22ll23ll . and can be easily extended from gamma processes to general completely random 
measures. 

We also proposed dynamical extensions of our model for both discrete and continuous time data, and 
applied it to modeling the bestsellers' lists on the New York Times. Our dynamic extension may be useful 
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Normalized weights 




Figure 3: Mean normalized weights for paperback nonfiction (left) and hardcover fiction (right). The 
black lines represent the weight associated to all the books that have not appear in the top-20 lists. 



for modeling time varying densities or clusterings as well. In our experiments we found that our model is 
insufficient to capture the empirical observation that bestsellers often start off high on the lists and tail off 
afterwards, since our model has continuous sample paths. We adjusted for this by simply not including 
books in the model prior to their publication date. It may be possible to model this better using models 
with discontinuous sample paths, for example, the Orstein-Uhlenbeck approach of l24l where the process 
evolves via a series of discrete jump events instead of continuously. 



Acknowledgements 

YWT thanks the Gatsby Charitable Foundation for generous funding. 



RRn° 8140 



Bayesian nonparametric models for ranked data 



12 



References 

[1] R.D. Luce. Individual choice behavior: A theoretical analysis. Wiley, 1959. 

[2] R. Plackett. The analysis of permutations. Applied Statistics, 24:193-202, 1975. 

[3] R.D. Luce. The choice axiom after twenty years. Journal of Mathematical Psychology, 15:215-233, 
1977. 

[4] D.R. Hunter. MM algorithms for generalized Bradley-Terry models. The Annals of Statistics, 
32:384-406, 2004. 

[5] I.C. Gormley and T.B. Murphy. Exploring voting blocs with the Irish electorate: a mixture modeling 
approach. Journal of the American Statistical Association, 103:1014-1027, 2008. 

[6] R Diaconis. Group representations in probability and statistics, IMS Lecture Notes, volume 11. 
Institute of Mathematical Statistics, 1988. 

[7] I.C. Gormley and T.B. Murphy. A grade of membership model for rank data. Bayesian Analysis, 
4:265-296,2009. 

[8] J. Guiver and E. Snelson. Bayesian inference for Plackett-Luce ranking models. In International 
Conference on Machine Learning, 2009. 

[9] F. Caron and A. Doucet. Efficient Bayesian inference for generalized Bradley-Terry models. Journal 
of Computational and Graphical Statistics, 21(1): 174-196, 2012. 

[10] G.R Patil and C. Taillie. Diversity as a concept and its implications for random communities. 
Bulletin of the International Statistical Institute, 47:497-515, 1977. 

[11] J. Pitman. Combinatorial stochastic processes. Ecole d'ete de Probabilites de Saint-Flour XXXII - 
2002, volume 1875 of Lecture Notes in Mathematics. Springer, 2006. 

[12] J. F. C. Kingman. Completely random measures. Pacific Journal of Mathematics, 21(l):59-78, 
1967. 

[13] M.K. Pitt and S.G Walker. Constructing stationary time series models using auxiliary variables 
with applications. Journal of the American Statistical Association, 100(470):554-564, 2005. 

[14] S. Watanabe. A limit theorem of branching processes and continuous state branching processes. 
Journal of Mathematics of Kyoto University, 8:141-167, 1968. 

[15] D. A. Dawson. Stochastic evolution equations and related measure processes. Journal of Multivari- 
ate Analysis, 5:1-52, 1975. 

[16] S.N. Ethier and RC Griffiths. The transition function of a measure-valued branching diffusion with 
immigration. Stochastic Processes. A Festschrift in Honour of Gopinath Kallianpur (S. Cambanis, 
J. Ghosh, RL Karandikar and PK Sen, eds.), 71:79, 1993. 

[17] R.H. Mena and S.G. Walker. On a construction of Markov models in continuous time. Metron- 
International Journal of Statistics, 67(3):303-323, 2009. 

[18] S.Feng. Poisson-Dirichlet Distribution and Related Topics. Springer, 2010. 

[19] J.C. Cox, J.E. Ingersoll Jr, and S.A. Ross. A theory of the term structure of interest rates. Econo- 
metrica: Journal of the Econometric Society, pages 385^1-07, 1985. 



RRn° 8140 



Bayesian nonparametric models for ranked data 



13 



[20] S. N. Ethier and R. C. Griffiths. The transition function of a measure-valued branching diffusion 
with immigration. Stochastic Processes, 1993. 

[21] L.F. James, A. Lijoi, and I. Priinster. Posterior analysis for normalized random measures with 
independent increments. Scandinavian Journal of Statistics, 36(l):76-97, 2009. 

[22] J.E. Griffin and S.G Walker. Posterior simulation of normalized random measure mixtures. Journal 
of Computational and Graphical Statistics, 20(l):241-259, 201 1. 

[23] S. Favaro and Y.W. Teh. MCMC for normalized random measure mixture models. Technical report, 
University of Turin, 2012. 

[24] J. E. Griffin. The Ornstein-Uhlenbeck Dirichlet process and other time-varying processes for 
Bayesian nonparametric inference. Journal of Statistical Planning and Inference, 141:3648-3664, 
2011. 



A Proof of Theorem 1 

The marginal probability ( fT~4| > is obtained by taking the expectation of ( fT3| l with respect to G. Note how- 
ever that ( fT3| is a density, so to be totally precise here we need to work with the probability of infinitesimal 
neighborhoods around the observations instead, which introduces significant notational complexity. To 
keep the notation simple, we will work with densities, leaving it to the careful reader to verify that the 
calculations indeed carry over to the case of probabilities. 



P{(Xt,Z t )ti) =E [P{{Y t ,Z t ) L l=l \G)\ 



=E 



A' 



e -c<x)£„z« l[G({X* k }) n *e 



G({X'})Y. H (S Uk -l)Zn 



k=l 



The gamma prior on G — Yl°jLi w j&Xj 15 equivalent to a Poisson process prior on TV = J2fLi ^(wj.Xj) 
defined over the space R + x X with mean intensity X(w)h(x). Then, 



=E 



K oo 

< ' - 1IX>'/' 1,A '< .v;;< • >: 

fe=U=i 



„- / wN(dw,dx) E« Zu 



ei(Stik—l)Zei 



Applying the Palm formula for Poisson processes to pull the k = 1 term out of the expectation, 



E 



-fw(N+8 w ,,)(dw,dx)T,u Z e , 



=E 



K oo 



K oo 

x {wl) nk h{Xt)e- w '^ 5ui -^ Zu X{w\)dwl 



„- / wN{dw,dx) E« z u JJ w ™« l(Xj = X* k )e- W ' ^ 

k=2j=l 

x h{Xl) ( {wl) nk e-<^ s ^ ZH \{wl)dwl 
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Now iteratively pull out terms k = 2, . . . , K using the same idea, and we get: 

K 



=E 



-G(X)E« z « 



k=l 



fe=i V « / 

This completes the proof of Theorem[T] 



(31) 



B Proof of Theorem 2 

Let / : X — > K be measurable with respect to ff. Then the characteristic functional of the posterior G is 
given by: 



E[e -//(x)G( rfa: )p ((y ^ Z£) L =i | G)] 

E[p((y,,^ =1 |G)] 



(32) 



The proof is essentially obtained by calculating the numerator and denominator of ( |32| >. The de- 
nominator is already given in Theorem [T] The numerator is obtained using the same technique with the 
inclusion of the term e-f f( x ) G ( dx ), which gives: 



E 

=E 



e-f f( - x)G{dx) P((Y e ,Z e )f =1 \G) 



K 



■/(/(*)+£« Z u )G(dx) 



k=l 



Uh(x* k ) /(^)^ e -:(/(^)+E 



Hw* k )dw* k 



By the Levy-Khintchine Theorem (using the fact that G has a Poisson process representation TV), 



; exp 



(1 - e 



-w(f(x)+J2 ei Zu) 



) X(w)h(x)dwdx 



K 

J] h(X* k ) / {wl)^e- w ^ x ^^ z ^\{wl)dwl 

k=l 



(33) 



Dividing the numerator ( pi) by the denominator ( |33"j ), the characteristic functional of the posterior G is: 



E 



cxp (- J(l - e- wf{x) )e~ E « z ^\{w)h{x)dwdx 



fe=i 



\(w* k )dw* k 



(34) 



Since the characteristic functional is the product of K + 1 terms, we see that the posterior G consists 
of K + 1 independent components, one corresponding to the first term above (G*), and the others cor- 
responding to the K terms in the product over k. Substituting the Levy measure X(w) for a gamma 
process, we note that the first term shows that G* is a gamma process with updated inverse scale t*. The 
kth term in the product shows that the corresponding component is an atom located at X k with density 
(w* k ) nk e~ Wk SeikZu X(wl); this is the density of the gamma distribution over w* k in Theorem]^ This 
completes the proof. 



RRn° 8140 



Bayesian nonparametric models for ranked data 



15 



C Proof of Proposition [4] 

We have 

P(G t {Xik) = Q\w t -i„k) = exp(-<l>t-iWt-i,k) 

Assume that 

P(G t (X lk ) = 0|w; sfc ) = exp(-y t \ s w sk ) 

then 



P(G t (Xik) = 0\w s -i ik ) = J ex-p(-y t \ s w sk )p(w sk \w s -i y k)dw sk 

= J2 /exp(- 2 / t|5 ^K fe |c s - 1 , fe )K Cs -i, fc h S -i, fe )^ 

Vt\s 



= ex p 



C,-l,k 



-C s _i,fcl0g 1 



p(c s -i >k \w s -i, k ) 



-yt\ s <Ps-i 

exp \ ; ; W s -l lk 

h-l + t + y t \ s 

D Gibbs sampler for the dynamic nonparametric Plackett-Luce model 

For ease of presentation, we assume that <fr t takes the same value cf> at each time step. The Gibbs sampler 
will iterate between the following steps 

1. a. Fort = 1, ... ,T, update G t (X) given (G t _i(X), a, (j>) 

b. Fort = 1, ...,T, update (c t ,c t *) given (w t ,w t *,w t+ i,w t+u , 0, a) 

2. a. Update a given (Z, <fi) 
b. Fort = 1,...,T 

Update given (c^i, , Z, 0, a) 
Update c t * given (w t *,Z, <j), a) 

3. Fort = 1, ... ,T, update (w t ,W t *) given (c t _i, c 4 _i*, c t , c t *, Z t , a, (j)) 

4. For t = 1, . . . , T, update Z t given (wt, Wt*) 

5. Update given w,w*,a 7 (f> 

The steps are now fully described. 
l.a) Sample (G t (X)) given (a, 0) 
We have 

Gi (X) | a ~ Gamma(a, r) 

and fort = 1,...,T- 1 

Gf+i (X) ~ Gamma(a + Mt , r + <j>) 
where M t ^Poisson(0Gt(X)). The weights (wt,wt*) are then appropriately rescaled. 
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l.b) Sample (c, c*) given (w, w*, <f>, a) 

Consider first the sampling of c\-t. We have, for t = 1, . . . , T and k = 1, . . . , K 

p(c t k\wtk,m+i,k) oc p(c tk \w t k)p(w t +iM c tk) 

where 

p{c t k\wtk) = Poisson(c tfe ;c/>w tfe ) 

and 



^ t+1>fc |c tfe ) = | Gamma (w t+ i 



if wtfe = 
fc ;c t fe,T + 0) if w tfe > 



mm 



Hence we can have the following MH update. If w t +\ t k > 0, then we necessarily have c t k > 0. We 
sample c* tk ^zPoisson(0ui t fe) where zPoisson(0u> t fc) denotes the zero-truncated Poisson distribution and 
accept c* k w.p. 

' Gamma(w t+ i,fc; c* tk ,r + 4>) ' 
' Gamma(w t+ i !fc ; c tk , r + <p) / 

If w t +i.k = 0, we only have two possible moves: c tfe = or = 1, given by the following 
probabilities 

P(c tfc = 0\w t+ i,k = 0, w 4 fe) = 



P(c tk = l\w t+ i,k = 0,wtk) 



exp(—<fw) t k) + 4>wtk exp(-0w tfe )(r + (j>) 1 + 0w*fc( r + 4 

(f>W tk exp(-(f)Wtk)(T + 4>) _ 4>W t k(.T + 4>) 



exp(-<jnv t k) + <t>w t k exp(-^w tfe )(r + </>) 1 + <jnv tk (T + 0) 



Note that the above Markov chain is not irreducible, as the probability is zero to go from a state 
(c t k > 0, Wt+i,k > 0) to a state (c t k — 0, w t +i,k = 0), even though the posterior probability of this event 
is non zero in the case item k does not appear after time t. We can resolve that by the following procedure, 
that uses a backward forward recursion. 

Assume that item k does not appear after time step r k (the same procedure applies if item k does 
not appear before time step t^T). Then we can sample jointly the whole sequence (wk,t, Ck,t)t=r k +i,...,T 
using the following backward forward recursion. 

Let 



x T = Z Tk (35) 

fe=i 



and fori = T- 1,...,t£ 



tk 

fc=l 



<t>X t +l 

- 4> + x t +i 



We have, for k = 1, . . . , K and t = t£ 



c tk \{Z,(t>, w tk ) ~ Poisson ( 1 ^ f — 4>w tk ) (36) 
V 1 + <p + x t J 

w t +i,k\ c tk, Z ~ Gamma (c M , r + <f) + x t+1 ) (37) 
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2.a) Sample a given (Z, 6) 

We can sample from the full conditional which is given by 

a\(Z,j,6) ~ Gamma (a + if, & + ?/! +log(l + zi)) (38) 
where x\ and yi are obtained with the following recursion 

m 

X T = Y. ZTk (39) 

k=l 

y T = (40) 



andfort = T- 



4>x t+ i 
- + 6 + x t+1 

1 + 6 



x t = Ztk + 7 
k=i ' 

yt = yt+i - log 



1 + 6 + x t+1 



2.b) Sample (c*, w*) given (Z, 0, a) 

We can sample from the full conditional which is given by 

wi* | (Z, 6, a) ~ Gamma (a, t + x\) (41) 
where x Y is defined above. Then for t = 2, . . . , T, let 

1 + 6 

— r-; <M-i* 



Ct-i*\(Z, 6,a,w t -u) ~ Poisson 

«>t* |c t -i* , a ~ Gamma (a + Ct_i* , r + + x t 



3) Sample (w,w*) given (Z, a, c, c*, 6) 

For each time step t = 1, . . . , T 

• For each item fc = 1, . . . ,K, sample 

w tfe |c t _i ife ,c tfe ,Z t - Gamma ^n tfe + Ct_i, fe + c tfe , r + 20 + ^ 5 tife Z ti ^ (42) 

if Ctfe + Ct-^fe + n t k > 0, otherwise, set w t k = 0. The occurence indicator S t ik is defined as 

(0 if3j<iwithF tj =X fe *; 
Otifc — < , . (43) 

1 otherwise. 
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Sample the total mass 

w t *\ct*,Ct-u,Z t ,a ~ Gamma a + c t * + Ct_i*, r + 2<j) + Z H (44) 



1=1 



4) Sample Z given (to, tu*) 

For i = 1, . . . , T and i = 1, . . . m, sample 



Z ti \w, — Exp + ^2 StikWtkj 



5) Sample 4> given tu, u?*, a, 



(45) 



We sample <j> using a MH step. Propose <j> = cj) exp(cre) where a > and e <~ A/"(0, 1). And accept it 
with probability 



P{wt+l*\<l>, w t*) TT g(Wt+^fc|0/Wtfc) 

p(tot+i*|0,iut*) ^ p(wt+i,fe|0,w tfe ) 



(46) 
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